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Abstract. We discuss an effective numerical method for 
solving the transport equations for cosmic ray nucleons in 
the atmosphere. It is demonstrated that the nucleon attenua- 
tion lengths are strongly energy and depth dependent due to 
the non-power-law behavior of the primary spectrum, growth 
of the total inelastic cross sections with energy, and scaling 
violation in the nucleon-nucleus interactions. The numerical 
results are compared with the available experimental data. 



1 Introduction 

Measurements of the fluxes of secondary cosmic -ray protons 
and neutrons can furnish valuable information about primary 
cosmic rays and about the nuclear interactions at high ener- 
gies. In order to extract this information from experimental 
data, it is necessary, among other things, to be able to calcu- 
late the nucleon attenuation lengths, which are functionals of 
the primary spectrum and which also depend, in general, on 
energy and atmospheric depth. 

In this paper we discuss some results obtained by using a 
simple but rather efficient and physically evident numerical 
method for solving transport equations describing the propa- 
gation of cosmic-ray protons and neutrons through the atmo- 
sphere (Naumov and Sinegovskaya, 2000). The method is 
applicable at sufficiently high energies and demands several 
approximations which are quite traditional for high-energy 
atmospheric cascade calculations. Namely, we use the one- 
dimensional (ID) approach, in which all secondary particles 
in the cascade are supposed collinear with the projectiles. 
The processes of generation of NN pairs in meson-nucleus 
collisions are disregarded, but the corresponding contribution 
is typically small and can be taken into account as a correc- 
tion (Vail et al., 1986). Besides, we apply the standard su- 
perposition model for the collisions of cosmic -ray nuclei and 
neglect the geomagnetic effects and proton ionization energy 
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losses. These approximations confine the range of applica- 
bility of the method. However, it remains rather broad in or- 
der to describe all currently available data on the high-energy 
nucleons in the atmosphere and at sea level. In particular, it 
covers completely the depth-energy area relevant to calcu- 
lations of atmospheric muon and neutrino fluxes at energies 
above several GeV, - one of the most important fields of ap- 
plication of the method. 

We compare our calculations with the experimental data 
and with the results of a more sophisticated approach by 
Fiorentini et al. (2001) based on an updated version of CORT 
code which was designed for applications to low and inter- 
mediate energy atmospheric cascade calculations and which 
takes into account essentially all significant effects (Naumov, 
1984; Bugaev and Naumov, 1985a,b; Bugaev et al., 1998). 



2 The Z factor method 

Within the above assumptions, the problem of calculating the 
differential energy spectra of protons D p (E, h) and neutrons 
D n (E, h) at the depth h consists in solving the following set 
of ID transport equations: 
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(N, N' = p, n), with the boundary conditions 

D p (E,0) = Dp(E), D n (E,0) = D° n (E), 



(1) 



(2) 



where D^{E) and D^(E) are the differential energy spec- 
tra of protons and neutrons at the top of the atmosphere; 
\ N (E) = l/[7Vocr^5}(£')] is the nucleon interaction length; 
d(JN'N(E', E)/dE is the differential cross section for inclu- 
sive reaction N'A — > NX; E' and E are the total energies 
of the projectile and final nucleons, respectively. 
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The approximate isotopic symmetry of NA interactions 
makes it possible to reduce the set of equations (1) to two 
independent equations for the linear combinations 

N±(E,h) = D p (E,h)±D n (E,h). (3) 

Let us write these combinations in the form 

N±(E,h) = N±(E, 0)cx P - — | - . (4) 

Below, the functions A±(E, h) will be referred to as effec- 
tive attenuation lengths. It is also convenient to introduce the 
dimensionless functions 



Z±{E,h) = l-X N {E)/k±{E,ti) 



(5) 



("Z factors"). In just the same way as the effective atten- 
uation lengths, the Z factors contain full information about 
the kinetics of nucleons in the atmosphere. Substituting Eqs. 
(3), (4), and (5) into transport equations (1) and integrating 
by part, we find that the Z factors satisfy the integral equation 

Z±{E,h) = \ [ dti [ dxr)±{x,E)®±(x,E) 
h Jo Jo 

xexp[-h'D±(x,E,ti)}, (6) 
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Although Eq. (6) is nonlinear, it is much more convenient to 
solve it by an iterative procedure than the original transport 
equations. 1 The rate of convergence of the procedure de- 
pends on the choice of zero-order approximation. The sim- 
plest choice is Z± ^ (E, h) = 0, in which case 
is independent of h and in first approximation we have 
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Obviously, the recurrent relations for n-th approximation (for 
n = 1, 2, . . . ) are given by 

Z ( £ ) (E,h) = \( dti [ dxr]±(x,E)<$>±(x,E) 
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Numerical analysis has revealed that the rate of convergence 
of this algorithm is quite sufficient for practical uses. 

'The main reason is in the comparatively weak energy depen- 
dence of the Z factors (see Fig. 1 below). 



3 Numerical results and discussion 

In the present calculations we employ three models for the 
primary cosmic ray spectrum and composition: NSU (Nikol- 
sky et al., 1984; Nikolsky, 1987), EKS (Erlykin et al., 1987), 
and FNV (Fiorentini et al., 2001). The NSU and EKS mod- 
els describe the high-energy range of the spectrum (~ 10 2 
to ~ 10 8 GeV/nucleon) and take into account the change 
of the spectral index in the "knee" region. The FNV model 
is a parametrization of the data from the most recent mea- 
surements of the primary spectrum below the "knee" (0.1 
to about 10 4 GeV/nucleon) for minimum of solar activity. 
All three spectra are extrapolated up to energy E = E c = 
3 x 10 10 GeV/nucleon above which a soft cutoff is supposed. 

For the differential cross sections d<TN<N(E', E)/dE we 
use slightly revised semiempirical model by Kimel' and Mo- 
khov (1974, 1975). For the total inelastic cross section a p n ^ 
we apply the parametrization by Mielke et al. (1994). 

Our calculations based on the Z factor method are per- 
formed for the energy range between 1 and 3 x 10 10 GeV at 
h < 4 x 10 3 g/cm 2 . It should be noted that, at h ~ 10 3 g/cm 2 , 
the energy losses by protons are important up to E ~ 30 
GeV. 2 The Z factor method (which neglects the energy loss 
effect) is extrapolated to the low-energy region in order to 
match the results of a more accurate analysis performed with 
CORT code. In the absence of experimental data on the 
fluxes of secondary nucleons arriving from oblique direc- 
tions, our calculations for h > 10 3 g/cm 2 are used at present 
only to test convergence of the iterative algorithm. 

At all values of E and h iterative process converges fast: 
five to six iterations are sufficient for calculating the Z fac- 
tors to precision not poorer than 10~ 3 — 10~ 4 . At moderate 
depths, h < 300 g/cm 2 , even the first approximation ensures 
a precision of a few percent, which is sufficient for many 
applications of the theory - in particular, for calculating the 
fluxes of atmospheric muons and neutrinos. 

Figure 1 illustrates energy dependence of the Z factors cal- 
culated using the EKS model of the primary spectrum for 
several oblique depths. The ground of the observed depen- 
dence of Z± on E and h is in the following three effects: (i) 
a non-power-law primary spectrum, (ii) energy dependence 
of the total inelastic cross section, (iii) violation of Feyn- 
man scaling in NA interactions. As a result, the attenuation 
lengths are also energy and depth dependent. This fact must 
be taken into account explicitly for determining the nucleon 
interaction length A^v (E) from the measured nucleon inten- 
sities in the atmosphere. Local minima in the Z factors that 
appear in the region around 45 GeV are due to the beginning 
of the growth of a^^E). At not overly large depths, the 
shape of the energy dependence visibly changes at E > 10 6 
GeV, which is caused by artificially introduced freezing of 
the growth of the quasielastic peak in the reaction pA — > pX. 
The vanishing of the Z factors at E — E c is due to the sup- 



2 By way of example, we indicate that, at E = 10 GeV, the rele- 
vant correction is about 20% for protons and about 6% for neutrons 
(Bugaev and Naumov, 1985a,b). 
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Fig. 1. Z+ and Z- vs energy, computed with the EKS primary 
spectrum for seven oblique atmospheric depths. Both Z+ and Z- 
decrease with increasing depth. 



posed cutoff of the primary spectrum at E > E c . 

Figure 2 shows the atmospheric growth of the proton fluxes 
for five momentum bins. The data of the balloon-borne ex- 
periment CAPRICE 94 (Francke et al., 1999) are compared 
with our calculations performed for the same bins by using 
the Z factor method and CORT code. In both calculations 
we use the FNV model of the primary spectrum. The nomi- 
nal geomagnetic cutoff rigidity in the experiment was about 
0.5 GV; hence the geomagnetic effect can be neglected. As 
one can expect, the proton fluxes calculated by the Z factor 
method systematically exceed those are obtained with CORT 
code, which takes into account the energy loss effect. The 
discrepancies become larger at increasing depth but vanishes 
at increasing the proton momentum. Minor differences be- 
tween two calculations at small depths are in part due to dif- 
ferent treatments of nucleus-nucleus interactions. 

Comparison of the calculated differential energy spectra 
of nucleons at three atmospheric depths with the data from 
Fawler et al. (1966); Aguirre et al. (1968); Apanasenko and 
Scherbakova (1968) are depicted in Fig. 3. The data col- 
lected in Fig. 3 were obtained indirectly, from measurements 
of photon spectra in extensive air showers (Grigorov et al., 
1973) and are therefore model-dependent to considerable ex- 
tent. 3 Nonetheless, our calculation relying on the EKS and 
NSU models of primary spectrum are by and large consistent 
with this data sample. 

In Fig. 4 we compare the calculated energy spectra of 
protons and neutrons at sea level with the data from Brooke 
and Wolfendale (1964); Ashton and Coats (1968); Ashton 
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3 In particular, approximate formulas used for recalculating pho- 
ton energy to nucleon energy lead to a distortion of the nucleon 
energy spectra. 
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Fig. 2. Atmospheric growth curves for protons. The data are 
from the CAPRICE 94 balloon-borne experiment (Francke et al., 
1999). The curves are calculated by CORT and by using the Z fac- 
tor method with the FNV model of primary spectrum. The legend 
shows the proton momentum bins and scale factors (in parentheses). 



et al. (1969, 1970); Diggory et al. (1974); Kornmayer et al. 
(1995). The results of Monte Carlo calculation by Lumme et 
al. (1984) is also shown. Direct measurements of the proton 
energy spectra at sea level are very fragmentary and one can 
speak only about qualitative agreement with the results of 
our calculations. Estimates reveal (see also Vail et al. (1986); 
Bugaev et al. (1998)) that the inclusion of processes of nu- 
cleon production in meson-nucleus interactions can increase 
the vertical flux of nucleons at sea level by no more than 10% 
at E = 1 TeV and by about 15% at E = 10 TeV, but this in- 
crease is much smaller, in either case, than the uncertainties 
in the NA cross sections and in the primary spectrum. Ex- 
perimental data on the neutron component at sea level are 
more detailed, but rather contradictory. The results of our 
calculations are in quite a good agreement with the data from 
recent measurements at prototype of the KASCADE facility 
in Karlsruhe (Kornmayer et al., 1995). As can be seen from 
Fig. 4, the neutron data by Kornmayer et al. (1995) below 
200-300 GeV are described by the calculation with the EKS 
and FNV models of primary spectrum somewhat better than 
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Fig. 3. Energy spectra of nucleons at three atmospheric depths. 
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Fig. 4. Energy spectra of protons and neutrons at sea level. 



by the calculation with the NSU model. 

It can be hoped that further experiments to study the nu- 
cleon component of secondary cosmic rays will allow a more 
detailed test of the method and of the models for the primary 
spectrum and for nucleon-nucleus interactions. 
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